Dynamic Stimulation of Quantum Coherence in Systems of Lattice Bosons 
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Thermal fluctuations tend to destroy long-range phase correlations. Consequently, bosons in a 
lattice will undergo a transition from a phase- coherent superfluid as the temperature rises. Contrary 
to common intuition, however, we show that non-equilibrium driving can be used to reverse this 
thermal decoherence. This is possible because the energy distribution at equilibrium is rarely optimal 
for the manifestation of a given quantum property. We demonstrate this in the Bose-Hubbard model 
by calculating the non-equilibrium spatial correlation function with periodic driving. We show 
that the non-equilibrium phase boundary between coherent and incoherent states at finite bath 
temperatures can be made qualitatively identical to the familiar zero-temperature phase diagram, 
and we discuss the experimental manifestation of this phenomenon in cold atoms. 
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A system of bosons confined to a lattice has long rep- 
resented an alluring opportunity to study the interplay 
between two phenomena at the heart of the academic 
and industrial interest in many-body quantum mechan- 
ics: particle tunneling and phase coherence. The Bose- 
Hubbard model (BHM) describing such systems in the 
tight-binding approximation is much richer than its sim- 
ple mathematical form betrays. It admits such novelties 
as dynamic localization [TJ [2] , photon-assisted tunneling 
[3l H], as well as an archetypal example of a quantum 
phase transition between superfluid and insulator-like 
states [5j [6]. Thermal fluctuations destroy these quan- 
tum effects, but it is comparatively unknown that delib- 
erately driving the system out of equilibrium can mod- 
erate or reverse entirely the destructive effect of raising 
the temperature. 

Despite its obscurity, it has been known since the 
1960 's that pushing a system out of equilibrium can en- 
hance its quantum properties. In 1966, Wyatt et. al. [7] 
showed that illuminating a microbridge could stimulate 
its superconductivity. Eliashberg explained this result in 
1970 by calculating the nonequilibrium quasiparticle dis- 
tribution induced by the radiation [8 . Blamire et. al. 
later demonstrated that superconducting transition tem- 
peratures could be enhanced in this manner by several 
times their equilibrium values [9] . More recently, the idea 
of non-equilibrium phase transitions has garnered inter- 
est in studies of optically trapped atoms [10] and induced 
topological structure [TT] . 

The reason for enhancement goes as follows. The quan- 
tum properties of a system depend on the energy distri- 
bution of excitations. It is easily verified that the equilib- 
rium distribution is rarely optimal for the enhancement 
of a chosen property. A brief survey of the model at equi- 
librium reveals this to be the case for the BHM. Indeed, 
an analogue of photon-assisted tunneling (PAT) [2JE] has 



already been theorized for the BHM at T = 0. However, 
a fully non-equilibrium treatment including the effects 
of temperature (and how they may be mitigated) is not 
known to us. 

In this Letter, we shall show that harmonically driving 
a system of lattice bosons connected to a thermal reser- 
voir can increase the region in parameter space where 
the quantum coherent phase exists. Even for finite tem- 
peratures of the bath, the phase diagram of the BHM 
can be made qualitatively identical to the T = di- 
agram. We shall demonstrate this by defining non- 
equilibrium correlation functions (a\(t)aj(t f )) within the 
Keldysh and Floquet formalisms [T2HT6] . Divergence of 
the real part of this quantity corresponds to infinite long- 
range correlations. This will define the phase boundary 
between superfluid and incoherent states for our non- 
equilibrium system. We shall find these functions per- 
turbatively [17j[T8] in the small quantity J/U and arrive 
at a Dyson equation that has both ordinary and entry- 
wise (Hadamard) matrix products. This novel structure 
will then be solved by column vectorization for the sta- 
tionary non-equilibrium correlation function. 

To see how superfluidity may be enhanced by means 
of a non-equilibrium pulse, let us briefly review the BHM 
at equilibrium. The hamiltonian of the BHM is 

Hq+Hj = y^ajaj (a\di - l) -\i a]aj- ^ Jjja\aj , 

i i ij 

(i) 

where nearest neighbor tunneling of strength J^- = J 
is assumed on a 2D square lattice with T <C U. This 
bath temperature models the coupling to an environ- 
ment [19 that dissipates energy. Let us first consider 
the T = case where fi/U is close to some integer M so 
that fi/U = M ± S for some < S < 1/2. Letting the 
tunneling J be infinitesimally small, the ground state is 
a Mott insulator. The energy gap for adding a parti- 
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cle or hole is proportional to S and (1 — S) respectively. 
As 5 is tuned to zero (unity), the state with an extra 
particle (hole) becomes degenerate with the Mott insu- 
lator. Thus, excess particles (holes) will be free to hop 
among sites with no energy barrier. At low tempera- 
tures, they will condense producing superfluidity [6]. As 
J is increased, the low lying excitations become long- 
range collective particle (hole) tunneling events between 
the system and reservoir. These events promote particle 
fluctuations, but they tend to stabilize the phase across 
sites through the number-phase conjugate relationship. 
The manifold where the energies needed for these long- 
range excitations vanishes defines the phase boundary [5] 
at T = 0. When T is finite, sites will have a thermal prob- 
ability p n — e~P £n of being occupied by n bosons. Sites 
will have different energies, and their phases will rotate 
at different rates. This aggravates the phase fluctuations 
that destroy superfluidity. The region in J vs. \i space 
that permits superfluidity is thus reduced (especially near 
integer values of ji/U) as the temperature is increased. 
We will see that the long-range phase-coherence demon- 
strated by our non-equilibrium correlation functions can 
be interpreted physically as the artificial closing of the 
Mott gap thereby allowing for the excitation of long- 
range tunneling modes. 

To effect an enhancement of superfluidity, we shall 
need a perturbation that pushes the system out of 
equilibrium. A thermal bath is also necessary both 
to counterbalance the heating induced by the pertur- 
bation as well as to ensure that the concept of tem- 
perature remains well-defined. We shall add three 
terms to the Hamiltonian so that H (t) = Ho + 
Hj + i^bath + #cou P + H v (t). The bath and cou- 
pling Hamiltonians H hath = J2ia £ ^ b la h ^ and H coup = 
a\a,i model the coupling of each site 
to an infinite bath of oscillator degrees of freedom such as 
the collective modes of a larger condensate in the which 
the lattice is immersed. Since the only crucial function 
of the bath is to balance the energy input from the driv- 
ing, many types of dissipation are possible and the re- 
sults of this paper are likely to extend to these alter- 
nate models. While future work is needed to verify this 
assertion, low frequency irradiation of a Josephson ar- 
ray should easily corroborate or deny this intuition ex- 
perimentally. We shall model the strength of our bath 
coupling [19 by a purely local and Ohmic parameter 
= rje a exp {— e a / A}. The driving term that will force 
a departure from equilibrium is given by 

Hy (t) = V a l a i cos (k ■ x» - (It), (2) 

i 

In practice, Hy (t) describes what is called a Bragg pulse 
[lOj [20j |21] . In the limit where energy differences between 
internal atomic levels are much larger than J7, J, and T, 
the potential in Eq. ([2| is created by the superposition 



of two lasers offset to each other in both frequency and 
wave vector. The spatial dependence of our perturbation 
is necessary for nontrivial results because a constant per- 
turbation simply multiplies the correlation function by a 
spatially constant phase factor [15 . This factor is irrele- 
vant to the tunneling of particles, and it can be gauged 
away by a time-dependent transformation of our opera- 
tors. The precise form of the spatial dependence in our 
system, however, seems not to be very important and 
other schemes are certainly possible [2 . We have chosen 
this one because of its experimental simplicity. 

It is difficult to say much about the nature of the inco- 
herent and coherent nonequilibrium phases or how they 
relate to the Kosterlitz-Thouless transition. Our focus 
will be to determine the boundary between the two non- 
equilibrium phases. To do this we shall approximate the 
correlation function (a\(t)aj(t f )) and look for divergences 
of this quantity over long distances. We shall do this 
perturbatively in the small quantity J/U following the 
general method in [T71 [18]. Anticipating a nonequilib- 
rium formalism and considering only the first order in 
the self-energy, the correlation function can be written 
as an infinite sum of simple chain diagrams defined on 
the foward-backward Keldysh contour, C. The evolu- 
tion on C is important in nonequilibrium problems be- 
cause it dispenses with the need to know the state of the 
system at t = oo for the calculation of expectation val- 
ues. The contour- time-ordering is accounted for by allow- 
ing green's functions to have a matrix structure [T2| [14] . 

That is, if we define G = ( ^ K ^ R J where A, K, R refer 

to advanced, Keldysh, and retarded Green's functions, 
then the non-equilibrium Dyson series can be written as 
a sum of matrix products of correlation functions. 

/CO 
dh $«/ (t,ti)G il ,-(ti,t'). 
-co 

(3) 

where gi 3 - refers to the correlation functions with re- 
spect to the hamiltonians Ho + i^bath + H coup + Hy (t) 
at bath temperature T. Because these hamiltoni- 
ans are just sums of single-site terms proportional to 
products of density operators rii = a\ai, they are 
easy to diagonalize in the occupation basis. Addi- 
tionally, the bath can be decoupled from the system 
[T4| [T9] via a canonical transformation ai — » e s die~ s = 
aie^^ S( 6 ;« _6iQ = aiXi that uses the time-derivatives 
of the transformed fields to cancel the coupling to the 
bath. Averages of the form (T c a\ (t) aj (£')) simply trans- 
form to (T c a\ (t) a 3 (t'))(T c X\ (t) X 3 (t')) = g %3 (t - t f ) o 
fij (t — t') where o denotes a Hadamard or entry wise 
product given by [A o Bj = A a pB a p where a,/3 des- 

ignate components in the 2x2 Keldysh space. The 
inclusion of the driving field Hy (t) produces a simple 
phase factor [15 equal to e <S[Bin(k.x 4 -nt)-Bin(k.x 4 -nt')] m 



3 



We conclude that the non-equilibrium function is 
a product of the time-dependent factors produced by 
#bath + H coup + Hy (t) and the bare function with re- 
spect to Hq. Expanding the non-equilibrium prefactor in 
terms of Bessel functions, the exact expression for g^ is 

9ij(t,t') = t re (f-t')o/(t-t') 



X 



CO 



mn= — oc 



v 



xx, (-0 ,[—*,]. (4) 

The functions / and g^ re are 2x2 matrices of the equi- 
librium correlation functions for the bath and the sys- 
tem given by Ho respectively. As they are equilibrium 
functions, they depend only on the difference of their 
time-arguments. All of the non-equilibrium information 
is stored in the expansion of the phase prefactor which 
depends on t and t' separately. This is the source of 
the Jo dependence of the tunneling renormalization 
familiar from studies of dynamic localization [2]. 

The non-equilibrium phase factor makes Gij a function 
of r = (t + t') 1 2 rather than simply of A = t — t f . How- 
ever, due to the time-periodicity of Eq. (pi), Gij (r, A) 
is a function of r only up to period 2ir/Q. This dis- 
crete time symmetry of the Hamiltonian allows us to de- 
compose every matrix function in Eq. ^ as Gij (t,t') = 
h Eiv etNQr I-oo e ~ l ^G l0 (a;, TV). Following the tech- 
nique illustrated in [15], equation ^ can now be writ- 
ten in terms of the functions g^ (a;, TV), but it will be 
burdensome to work with it because the equation for 
Gij (cj, TV) will include contributions from Gij (c^, TV') for 
all TV' . Fortunately, we can mathematically represent 
this coupling as simple matrix multiplication if we trans- 
form to the so-called Floquet representation [15] [16] de- 
fined as Gij (cj) mn = Gij (uu + -^-^1), m — n). We may 
now think of Gij (cj) mn as an infinite square matrix in 
the two-dimensional space of Floquet indices m and n. 
Each element of this matrix is itself a 2 x 2 Keldysh ma- 
trix. We will suppress the site indices and make use of 
the discrete translational symmetry of the problem by 
transforming to lattice-momentum space G^(a;,q) mn = 
e l€li '( Xj ~ Xi> )Gij (<u) mn - Finally, we explicitly account 
for the plane-wave contribution to the full Green's func- 
tion by defining e^ - ™^"*^ (cj, q) mn = G;(cj,q) mn . 
Having made these transformations, we arrive at an ex- 
tremely simple form for the nonequilibrium Dyson equa- 
tion. 



0(w,q)=£M 1 - J(q,k) oQ(u,q) 



(5) 



where for given matrices A and B in the Floquet space 
(indexed by m,n), AB indicates an ordinary matrix 
product while A o B denotes a Hadamard in the Floquet 
space rather than Keldysh space. The matrix J (q, k) 



is a generalized lattice dispersion in Floquet space given 
by Jmn (q) = JY,v cos + {m-n) k v ] where v = 1, 2 
denotes principle directions in our square lattice. 

The existence of a Hadamard product in Eq. (|5| com- 
plicates matters. We cannot merely multiply by inverses 
to solve for Q because there are now two types of in- 
verses corresponding to the two types of products. This 
double-product structure in a Dyson equation seems so 
far unknown in any other context, but the situation can 
be salvaged by column vectorization (CV): the mapping 
of matrix A to a vector A consisting of the first column 
of A stacked on the next column and so on. We can then 
make use of the convenient identities relating Hadamard 
and ordinary matrix products through CV to solve Eq. 
([5| and rewrite it as 



g(u,q) = (l n >< - [l«p> 



>$(w)]£>{j(q)}) 1 g(oj), 
(6) 



where (8) indicates a Kronecker product, and D |^| de- 
notes the diagonal matrix with entries given by those of 
A. The identity matrix of size k x k is given by l kxk ^ 
while n p is the size of the matrix g in the Floqet (m, n) 
space. It signifies how many higher harmonics we wish 
to include, or equivalently, the maximum time-resolution 
of our treatment. If we needed infinite time-resolution, 
we would of course let n p — >• oo. However, each off- 
diagonal element (m, n) will be weighted by (J/U) m ~ n 
while g mn — » with increasing m + n, so we expect that 
rip need not be large to capture the relevant stationary 
behavior. 

To determine the phase boundary, we are only inter- 
ested in the stationary behavior of the system given by 
Qoo. Inverting the block-diagonal x matrix in Eq. 
(J6|, the real part of the Keldysh component ofour corre- 



lation function, Re^oo, can be displayed. Fig. [1(a)] shows 



the system at equilibrium (V = 0) including the effects 
of Ohmic dissipation. The phase boundary is given by 
the points where Re^oo diverges, and our results match 



those of Ref. [193 . Fig. |l(b) is an example of dynamic 
enhancement of superfluidity. Note the similarity of Fig. 
1(b) | to the T = phase diagram at equilibrium. In Fig. 
l(c)| we see the effect of making the pulse energy TiVi ten 
times smaller. Again there is supernuid enhancement, 
but notice that the valleys come to much finer points 
closer to the ji/U axis. This advocates an interpreta- 
tion wherein our perturbation excites phase-stabilizing 
collective modes and artificially closes the energy gap. A 
smaller perturbation energy TiVt is resonant with a smaller 
gap. Thus, the phase boundary is much closer to the fi/U 
axis at integer values of fi/U where the gap goes to zero. 
It becomes qualitatively identical to the T = equilib- 
rium phase diagram. 

The system's behavior in other regions of V, O space 
are catalogued in Fig. [2] For instance, at high frequencies 
(HQ ^> J7), we see the familiar phenomenon of dynamic 
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FIG. 1: Numerical density plots of Re^^ for bath temperat ure 
(3U = 25, coupling width A/U = 3, and strength rj = 0.01. [(a)] 
Equilibrium: V = 0. |(b)| Superflui d en hancement: V/U = 0.1, 



0.05, k = ^x, n p = 
V/U = 0.1, m/U = 0.005, k = 
the T = equilibrium diagram. 



(c) Superfluid enhancement: 
= 5. Note the similarity to 



localization [TJ [2]. All Wigner components other than 
= can be ignored, and the tunneling is renormalized 
J —> J Jo (V/£l) by the zeroth Bessel function of the first 
kind coming from the expansion of the non-equilibrium 
phase factor in Eq. Q. If V/Vt is tuned to a zero of 
Jo (x), we have dynamic suppression of tunneling. We 
find a featureless phase diagram (not shown in Fig. [2| 
because there is no region of J vs. /i space that admits 
long-range phase-coherence. 

This phenomenon, familiar from driven Josephson ar- 
rays, can be understood as a cancellation of the dynamic 
phase acquired due to one period of driving with that 
of the hopping between sites. The bosons become local- 
ized with no long-range correlations, and Re(?oo diverges 
nowhere. If we now tune V toward zero, the phase bound- 
ary returns to its equilibrium form. If instead Vt is tuned 
lower, it will eventually be small enough to be resonant 
with the low energy modes available when fi/U is close 
to an integer. We will again have dynamic enhancement 
of superfluidity. However, further from the valley, there 
will be no available low-energy modes, and we will see 
suppression of tunneling. The result is larger Mott lobes 
that almost touch the ji/U axis. 

We have demonstrated the enhancement of the super- 
fluid region in parameter space by driving. The exper- 
imental signature of this effect is similar to what has 
been found in time-of-flight experiments [22]. When \i 
and J are tuned to a point within the enhanced super- 
fluid region close to integer values of fi/U, there will be 
well-defined peaks in momentum space when the pertur- 
bation is on (V 7^ 0) and a featureless interference pat- 
tern corresponding to destroyed phase coherence when 
the perturbation is off (V = 0). 
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FIG. 2: The T = 0.04*7 equilibrium phase boundary (dotted line) 
plotted with examples from different regions of driving parameter 
space. 
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